Visualizing Spatial Data in R

library(rgdal)
# streams <- readOGR(dsn = "../data/R_GIS_Layers/", 
#                    layer = "Cove_Drainage_WGS84")
summary(streams)

plot(streams, col="blue")

## note that plots using base R may be distorted.

library(ggplot2)
streams_gg <- fortify(streams)  ## or use broom::tidy()
# str(streams_gg)
ggplot() + 
  geom_path(data = streams_gg, 
            aes(x = long, y = lat, group = group), 
                color = "blue", size = 1.5)

## note that without a coordinates correction, the plot may be distorted.

library(ggmap)
# Might need older version of ggplot2:
# devtools::install_github("hadley/ggplot2@v2.2.0")
# library(ggplot2)
myLocation <- c(lon = -109.2279, lat = 36.545) 
myMap <- get_map(location = myLocation, source = "google", 
                 maptype = "satellite", crop = FALSE, zoom = 13)
g <- ggmap(myMap, darken = c(0.25, "white"))
g + 
  geom_path(data = streams_gg, 
            aes(x = long, y = lat, group = group), 
            color = "lightblue", size = 1.5)
ggsave("shp_and_ggmap.png", dpi = 72) 

library(dplyr)
alum_data <- plot_data %>% filter(Analyte == "ALUMINUM")
g + 
  geom_path(data = streams_gg, 
            aes(x = long, y = lat, group = group), 
            color = "lightblue", size = 1.5) +
  geom_point(data=alum_data, 
             aes(x = Longitude, y = Latitude, size = Result, fill = Result),
             shape = 21, alpha = 0.9) +
  scale_fill_gradient(low = "#1f77b4", high = "#d62728") +
  scale_size(range = c(2, 8)) +
  guides(size = FALSE) +
  coord_equal() 
ggsave("static_map.png", dpi = 72) 

head(alum_data)

library(scatterplot3d)
with(alum_data, {
  lollipop_plot <- scatterplot3d(streams_gg$long, streams_gg$lat, 
                                 rep(0, length(streams_gg$long)),  # x y and z axis
                              color = "blue", pch = 16, type = "p",
                              angle = 120,            
                              scale.y = 1.75,           
                              main = "Example Lollipop plot",
                              xlab = "Longitude",
                              ylab = "Latitude",
                              zlab = "Concentration values",
                              xlim = c(-109.273, -109.15),
                              ylim = c(36.50, 36.59),
                              zlim = range(0, 46000),
                              grid = TRUE,
                              box = FALSE)
  
  # add the legend
  legend("topright", inset = 0.07,      
         bty = "n", 
         title = "Type of Results",
         c("High","Low"), 
         fill = c("#1B9E77", "#D95F02"))
  
  # add the lollipop points
  lollipop_plot$points3d(Longitude, Latitude, Result,        
                      col = "#282830", 
                      pch = 21, 
                      bg = ifelse(Result > 10000, "#1B9E77", "#D95F02"), 
                      lwd = 1,        
                      type = "h",      
                      cex = (Result / 50000) + 1)
})

library(rgl)
library(dplyr)
plot_3d <- with(alum_data, 
  plot3d(Longitude, Latitude, log(Result),        # x y and z axis
         col = ifelse(Result > 10000, "#1B9E77", "#D95F02"), size = 5,
         type = "p"))
rglwidget(elementId = "plot3drgl") # to show in presentation

## plotly recommends the development version of ggplot2, but will also work with the
## latest version of the package, version 2.2.1 
# library(devtools)
devtools::install_github('hadley/ggplot2')
library(ggplot2)
library(plotly)
p <- plot_ly(alum_data, 
             x = ~Longitude, y = ~Latitude, z = ~Result,
             marker = list(color = ~Result, 
                      colorscale = c('#FFE1A1', '#683531'), showscale = TRUE)) %>%
  add_markers() 

p

interactive_map <- g +  
  geom_point(data = alum_data, 
             aes(x = Longitude, y = Latitude, size = Result, fill = Result),
             color = "#000000", shape = 21, alpha = 0.8) +
  scale_fill_gradient(low = "#1f77b4", high = "#d62728",
                      guide = "legend") +
  scale_size(range = c(2, 8)) +
  coord_equal() +
  ggtitle("ALUMINUM")

ggplotly(interactive_map + theme_void(), filename="plotly")

library(animation)
ani.record(reset = TRUE)
for(a in unique(plot_data$Analyte)) { 
  plot_data_a <- subset(plot_data, Analyte == a)
  animated_maps <- g + 
    geom_path(data = streams_gg, 
              aes(x = long, y = lat, group = group),
              size = 1.5, color = 'lightblue') +
    geom_point(data = plot_data_a, 
               aes(x = Longitude, y=Latitude, size = Result, fill = Result),
               shape = 21, alpha = 0.8) +
    scale_fill_gradient(low = "#1f77b4", high = "#d62728") +
    scale_size(range = c(2, 8)) +
    guides(size = FALSE) +
    coord_equal() +
    ggtitle(a)
  
  print(animated_maps)
  ani.record()
}
oopts = ani.options(interval = 1)
ani.replay()
saveHTML(ani.replay(), img.name = "animation_plot", 
         htmlfile = "animation_5metals.html", ani.width=600, ani.height=400)

# Show neighbors on map of subbasins
library(sp)
class(lulc.sp1)

spplot(lulc.sp1, zcol = "PCTHERB", col = "black", 
       main = "Percent Herb")

library(spdep)

## determine who shares a border
ctchmt.nb1 <- poly2nb(lulc.sp1, queen = TRUE)
# put neighbors into list
ctchmt.nbwts.list <- nb2listw(ctchmt.nb1, style = "B")

moran.test(lulc.sp1$PCTHERB, ctchmt.nbwts.list)

par(mar=c(5, 5, 2, 2))
moran.plot(lulc.sp1$PCTHERB, listw = ctchmt.nbwts.list, 
           labels = lulc.df1$PCTHERB,
           xlab = "PCTHERB", ylab = "Spatially Lagged PCTHERB",
           main = "Moran Plot for PCT HERB")

## The network graph may require the following graph package that need to be 
## installed from source.  (This package is not generally required to create  
## similar graphs; they are only required for the example.  See the tutorial for  
## other examples based on csv files.)
# source("https://bioconductor.org/biocLite.R")
# biocLite("graph")
# biocLite("Rgraphviz")

library(magrittr)  ## for %<>% function
library(visNetwork)  ## to create the interactive network graph
library(jsonlite)
library(igraph)
network_interactive <- function(g, nodes, coords){
  ## first, switch graph construction to data frames:
  nNodes <- length(nodes)
  nodes_df <- matrix(NA, nrow=nNodes, 4)  # save first 4 attributes
  for (i in 1:nNodes){
    nodes_df[i, ] <- unlist(jsonlite::fromJSON(nodes[[i]]))[1:4]
  }
  nodes_df <- as.data.frame(nodes_df, stringsAsFactors=FALSE)
  
  ## because edges uses the names, need id to equal names of nodes, not numbers
  names(nodes_df)[1:3] <- c("old_id", "id", "type")
  
  ## and id column needs to go first: (only for igraph package, not visNetwork)
  ## igraph packge was used to help decide node coordinates
  nodes_df <- nodes_df[, c('id', 'old_id', 'type')]

  g_edges <- Rgraphviz::buildEdgeList(g)
  nEdges <- length(g_edges)
  edges_df <- NULL
  for (i in 1:nEdges){
    tmp <- c(from=g_edges[[i]]@from,
             to=g_edges[[i]]@to,
             unlist(g_edges[[i]]@attrs))
    edges_df %<>% bind_rows(tmp)
  }
  edges_df$weight <- as.numeric(edges_df$weight)
  
  ## Add layer info/ color column
  nodes_df$col_layer <- 
    ifelse(nodes_df$id %in%  grep("MainInput", nodes_df$id, value=T), 1, # 
           ifelse(nodes_df$id %in% grep("OtherInput", nodes_df$id, value=T), 3,
           ifelse(nodes_df$id %in%  grep("Midvalue", nodes_df$id, value=T),  4,
           ifelse(nodes_df$id %in%  grep("MidEqn", nodes_df$id, value=T),  5,
                  6))))
  
  ##  add coordinates to determine layout of each node:
  nodes_df <- full_join(nodes_df, coords)

  ## specify colors to use for each col_layer
  graph_colors <- c("gold", "darkorange", "tomato", 
                    "palegreen", "seagreen", "royalblue")
  # frame = borders
  frame_colors <- c("gold3", "darkorange3", "tomato4", 
                    "palegreen3", "seagreen4", "royalblue4")
  ## (visNetwork doesn't like color names with numbers)
  frame_colors_rgb <- rgb(t(col2rgb(frame_colors)), maxColorValue = 255)

  ## Add attributes so that the graph looks good:
  nodes_df$shape <- "ellipse"
  nodes_df$color.background <- graph_colors[nodes_df$col_layer] 
  nodes_df$color.border <- frame_colors_rgb[nodes_df$col_layer]
  nodes_df$color.highlight.border <- "darkred"
  
  edges_df$arrows <- "to"  # draw arrows on the edges.
  edges_df$color.highlight <- "black"
  
  network <- visNetwork::visNetwork(nodes_df, edges_df, width="100%", physics = TRUE) %>%
    visNetwork::visIgraphLayout() %>%
    visNetwork::visOptions(highlightNearest = TRUE, nodesIdSelection = TRUE,
               manipulation = FALSE) %>%
    visNetwork::visEdges(color = "slategray", smooth = list(roundess = 1)) %>%
    visNetwork::visNodes(font = list(size = 22), shape="ellipse")

  return(network)
}

network_interactive(g = networkDag, nodes = networkNodes, coords = networkCoords)

par(mar=rep(0, 4))  ## change margins so plot fills entire space
## default colors:
plot(1:8, pch=16, cex=5, col=1:8)

## Defining a new palette
new_palette <- c("darkred", "chartreuse", "turquoise", "purple",  
                 "gray45", "plum", "black", "#F08800")
palette(new_palette)
plot(1:8, pch=16, cex=5, col=1:8)

palette("default")  ## return palette to default colors

sessionInfo()